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Abstract. Stochastic models of chemical systems are often analysed by solving the correspond- 
ing Fokker-Planck equation which is a drift-diffusion partial differential equation for the probability 
distribution function. Efficient numerical solution of the Fokker-Planck equation requires adaptive 
mesh refinements. In this paper, we present a mesh refinement approach which makes use of a 
stochastic simulation of the underlying chemical system. By observing the stochastic trajectory for 
a relatively short amount of time, the areas of the state space with non-negligible probability density 
are identified. By refining the finite element mesh in these areas, and coarsening elsewhere, a suitable 
mesh is constructed and used for the computation of the probability density. 



1. Introduction. Stochastic simulation algorithms (SSAs) have been success- 



fully used in recent years to understand a number of biochemical models [3}[3Q, 39 
However, a systematic analysis of these models is challenging because of the compu- 
tational intensity of SSAs. A suitable alternative to stochastic simulation is a solution 
of the chemical Fokker-Planck equation 1 1^|25] . Consider a well-mixed system of N 
chemical species and denote by x = {xi^ . . . ,xn) a vector of concentrations of these 
species. The stationary chemical Fokker-Planck equation is a drift-diffusion partial 
differential equation (PDE) for an A^-dimensional probability distribution function 
p = p(x.) where x G 1^ C M^, which can be written in the following form: 

div [D(x)Vp(x) - v(x)p(x)] = (1.1) 

where D = D{x.) : Q ^ R^""^ is the diffusion matrix, v = v(x) : 1] ^ is the drift 
term and div is the divergence operator. There have been several methods developed 
in the literature to solve ( |1.1| ) for moderately large N. They include adaptive finite 
element methods (FEMs) which are commonly used for A/" < 3 4 and sparse grid 
approaches which are applicable for larger values of N 40 . Adaptive FEMs can be 
used to identify a suitable mesh which is refined in crucial regions and not in others. 
Although these approaches can be used for the chemical Fokker-Planck equation, they 
do not exploit the fact that its solution is a probability distribution of a stochastic 
process. 

In this paper, we present an adaptive mesh construction which is suitable for the 
FEM solution of ( |1.1| ) if this equation arises from modelling of stochastic chemical 
systems. The main idea is to exploit the fact that stochastic trajectories spend a 
significant amount of time in parts of the state space where the mesh refinement is 
needed. Since adaptive FEMs are mostly applicable for systems up to = 3, we will 
focus on systems of 3 chemical species. However, the presented methodology can be 
modified for larger (multiscale) chemical systems provided that they have up to n < 3 



important (slow) variables. In 13 20j, a method for estimation of coefficients of an 
effective Fokker-Planck equation is presented. This effective equation is of the same 
form as ( |1.1[ ) but it is written in the dimension n which is smaller than a total number 
A^ of chemical species. If one has a suitable SSA for simulating the low dimensional 
slow dynamics [9, 10 17 , the presented mesh refinement can be applied. However, 
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[1] Calculate propensity functions a/c(X(t)), /c = 1, 2, . . . , M. 



[2] Waiting time r till next reaction is given by (2.1). 

[3] Choose one j G {1,2,..., M}, with probability aj(X(t))/ao(X(t)), and perform 
reaction Rj^ by adding Vji to each Xi{t) for alH = 1, 2, . . . , A^. 

[4] Continue with step [1] with time t = t-\-r. 

Table 2.1 
The pseudo code for the Gillespie SSA. 



since the main aim of this paper is to present this novel numerical methodology, we 
will not consider any dimensional reduction and restrict to systems which are directly 
written in terms of N = 3 chemical species. In the following paper 12 , we will show 
how this mesh refinement can be used to study bifurcation behaviour of stochastic 
chemical systems. 

The paper is organised as follows. In Section [2j we introduce our notation, the Gille- 
spie SSA and the chemical Fokker-Planck equation. In Section [3j we introduce the 
standard finite element framework that we will use throughout the paper. We also 
formulate the problem in its weak form in order to define the problem to be solved. 
In Section |4j we introduce the stochastic simulation assisted adaptive Finite Element 
Method (saFEM). In Section [sj we offer some insight into how the algorithmic param- 
eters of the saFEM may be chosen in practise. In Section [6j implementational issues 
related to the algorithm are addressed. Numerical results concerning convergence of 
the method are presented in Section [71 

2. Chemical Fokker-Planck Equation. Let us consider a well- mixed system 
of N chemical species X^, i = 1,2, ...,A/', which are subject to M chemical reac- 
tions Rk, k = 1,2,..., M. The state of the system is described by the state vector 
X(t) = X2(t), . . . ,XAr(t)] where Xi{t) denotes the number of molecules of the 

corresponding chemical species. Each chemical reaction is described by its propensity 
function and the stoichiometric vector [IQ^,^^^- The propensity function a/c(x) is de- 
fined in such a way that afe(x)dt is the probability that the k-th reaction occurs in 
the infinitesimally small interval [t,t-\-dt) provided that X(t) = x. The stoichiometric 
vector is Uk = [z^/ci,^/c2, • • • ,z^/cAr] where Uki is the change in Xi during reaction Rk. 
Stoichiometric vectors form the corresponding stoichiometric matrix u = {i^ki)^i=i- 
The time evolution of the state vector X(t) is often simulated by the Gillespie SSA |23] 



which is described in Table 2.1 Given the values of the propensity functions (step 



[1]), the waiting time to the next reaction is given by: 

^^-^^T^V ^^^^^ ao(X(t)) = ^a,(X(t)), (2.1) 
<^o(X(t)) ^ 

and 7x is a uniformly distributed random number in (0, 1). The reaction Rj is chosen 
in step [3] using another uniformly distributed random number. 

The stationary probability distribution corresponding to the chemical system can be 
approximated by solving the chemical Fokker-Planck equation 25 . This equation can 
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be written in the form (1.1) where the diffusion and drift coefficients are: 



1 ^ 

dij{x.) = -^Vki^kjOtk{y^), z,j = 1,2, ...,Ar, (2.2) 



/c=l 

M N 



^,(x) = ^k^Oik{^) - ^ i = 1, 2, . . . , iV. (2.3) 

This equation is solved on a bounded domain Vt C [0, oo)^ in which the vast majority 
of the invariant probabihty density sits 18 . On the boundary dVt we introduce the 
homogeneous Neumann boundary condition: 

[D(x) Vp(x) + v(x)p(x)] • n(x) = 0, for x G dVt, (2.4) 

where n(x) is the outward facing normal at x G dO.. We seek a solution of 
which corresponds to the probability distribution function. Therefore, we impose the 
following normalisation condition: 

/ p(x)dx = 1. (2.5) 

The choice of the computational domain ^ might be problematic if we have no a 
priori information about the problem with (2.4) and (2.5). In this case, the 

stochastic simulations provide a reliable tool to determine the correct Vt as we will 
show in Section |4.2[ However, before discussing the details of saFEM, we have to 
introduce some finite element terminology. This will be done in the next section. 

3. Finite Element Method. Equation with boundary condition (2.4) can 
be numerically solved by the FEM. Since the finite element formulation is based on 
the corresponding weak formulation, we first introduce the weak solution p G H^(yt) 
by the equality 

a(p,(^) = 0, y^eH\n). (3.1) 

The bilinear form a(-, •) is naturally given by 

a(p, (/))=/ (I^(x) Vp(x) + p(x) v(x)) • V(/)(x) dx. (3.2) 



The finite element formulation is obtained by projecting the weak formulation (3.1) 
into a finite dimensional subspace Vh of H^{^}). Thus, we seek ph G Vh such that 

a{ph,cl)h)=0, ycl^heVh. (3.3) 

We note that taking a basis 0i, ^2, • • • , of Vh, we can express the finite element 
solution as 



Ph{x.) = ^pVi(x), X G r2. 



i=l 

Here, the coefficients G M solve the system of linear algebraic equations 



Ap = 0, 
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(3.4) 



where p = (jp^ -tP^ • • • ^p^)^ and the stiffness matrix A G M^><^ is defined by its 
entries 

Aij = a{(t)j,(t)i), ij = 1,2, ...,m. (3.5) 

We construct the finite element space Vh and the corresponding basis functions in the 
standard way 11 . In what fohows, we wih focus on computations in three dimensions, 
i.e. N = 3. We consider the finite element mesh Th consisting of tetrahedral elements. 
The lowest-order finite element space Vh then consists of globally continuous and 
piecewise linear functions over the mesh Th- 

Vh = {<Ph e H\n) : <l>h\K e F\K) for all K € %}, (3.6) 

where P^(i^) stands for the space of linear functions over the tetrahedron K G Th- If 
G M^, j = 1, 2, . . . , m, stand for the nodes of the mesh Th then the standard finite 
element basis functions are uniquely determined by the condition 

M^j) = ^ij^ hj = l,2,...,m, 
where 5ij stands for Kronecker's symbol. 



An efficient solution to problem (3.3) can be obtained by employing adaptively refined 
meshes 38 . The optimally adapted mesh leads to an approximation with the smallest 
error, provided the number of degrees of freedom is fixed. Practically, the optimal 
mesh can be hard to determine, but meshes close to the optimal can be found. These 
meshes are fine in regions where the solution exhibits steep gradients, boundary layers, 
interior layers or singularities, and they are relatively coarse in the other regions. 
The standard numerical approach for construction of nearly optimal meshes is the 
adaptive algorithm based on suitable a posteriori error estimators |4)[38]. This algo- 



rithm starts with an initial coarse mesh and refines it adaptively by a sequence of 



refinement steps. In each step, problem (3.3) has to be solved on the actual mesh, an 



error indicator has to be computed for each element, and based on these indicators the 
mesh is refined at suitable places. Using the mesh refinements assisted by stochastic 
simulations, we can avoid the sequence of refinement steps and construct a suitably 
adapted mesh at once. 

4. Adaptive mesh refinement assisted by stochastic simulations. Since 
the chemical Fokker-Planck equation is a continuous approximation of the evo- 
lution of the probability density given by the Gillespie SSA 23 , one can expect that 
trajectories simulated from the SSA will be informative. Once the trajectory has 
reached probabilistic equilibriunj^ regions surrounding the path of the trajectory are 
likely to be regions with non-negligible invariant density with respect to the steady 
state Fokker-Planck equation. We should be aiming to refine the finite element mesh 
in regions where the rate of change of the gradient of the invariant density is larger. 
The area in which we should be trying to refine our mesh can be well approximated 
by the region which has non-negligible invariant density. Therefore, stochastic simula- 
tions of the chemical system can be informative about a good choice of finite element 
mesh. 

The construction of the locally adapted mesh assisted by stochastic simulations is done 
in three stages. In Stage I, we use the stochastic simulations to identify those regions 



■"^Equilibrium can be reached simply by running the SSA until the state of the trajectory is 
sufficiently decorrelated from its starting position, or by starting at a steady state of the mean field 



approximation (4.1) of the chemical system. 
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[1] Identify steady states (stable and unstable) of the mean field approximation 
of the system, or in the case of oscillatory systems, one (or more) coordinates 
along the limit cycle. We denote these 5* G N points by {yk}k=i' 

[2] Run one (or more) Gillespie simulations of length B-\-T > for each of the initial 
conditions {yk}k=i' Here, B > is the length of a possible initial transient 
and T > 0. Take a subsample of Q > points per unit time per trajectory 



in the time interval [B^B + T] and denote them as in (4.2). We also denote 
by xf^^^ (resp. xf^^^) the maximal (resp. minimal) value of the i-th chemical 
species, i = 1, 2, 3, during time intervals [B^ B -\-T] of all S simulations. 

[3] Use (4.3), to define a neighbourhood F C [0,oo)^ of the points (4.2) as the 



region of the domain in which we require the finest level of refinement of the 



mesh. This is made up of ellipsoids around each point with radii given by (4.5) 
with parameter /3i > 0. 



[4] Define the domain of solution Q by (4.6), using parameter P2 > 



[5] Start with the mesh as one single cuboid covering the whole of Q. 

[6] Loop over all cuboids in the mesh. If the cuboid is sufficiently close (as given by 
( |4.7| )) to r, then split the cuboid into eight equally sized cuboids. Update the 
list of hanging nodes/hang type (face/edge hanging node as shown in Figure 
lOKleft)). 

[7] Repeat step [6] until the maximum resolution has been reached after H G N 
iterations, or the maximum total number of cuboids has been reached. 

Table 4.1 

Mesh generation using SSA trajectories. 



of the state space, where the system spends most of its time, and hence where we are 
interested to resolve the problem with the highest accuracy. In Stage II, we identify 
the computational domain as a cuboid that covers the region from Stage I and 
its neighbourhood. In Stage III, we construct the actual mesh in the computational 
domain Q based on the information from Stage I. In what follows, we provide detailed 
description of these three stages in the case of N = 3 chemical species. Following the 
notation of Section |2j the chemical species will be denoted as Xi, X2 and X3. The 
mesh generation algorithm is summarised in Table \AA] 



4.1. Stage I: Stochastic simulation. In step [1] of Table 4.1, we identify 
structures that a priori should exist in the probability density. An indication of the 
regions which may contain significant amounts of invariant density can be found by 
analysing the 3-dimensional system of ordinary differential equations (ODEs) 

= Vi{xi,X2,xs), {1,2,3}, (4.1) 

which is closely related the mean field approximation of the chemical system. Steady 



states of ODE system (4.1 ) will often coincide with regions of the solution of the steady 
state Fokker-Planck equation ( |1.1[ ) that have large density. The most important things 
to identify are stable steady states of ( |4.1[ ), which can be found by solving an algebraic 
system ['^i(xi, X2, X3), '^2(^1, ^2, ^3)7 '^3(^1, ^2, ^3)] = 0- In oscillating systems, one can 
identify limit cycles. There are several tools in the literature for analysis of ODEs of 
the form (|4l]), such as AUTO [15). 
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In step [2], these steady states can then be used as starting points for SSA trajecto- 
ries, once the numbers of molecules of each species have been rounded to the nearest 
integer. One might additionally wish to ensure that each chain is starting in proba- 
bilistic equilibrium by running a short simulation of length B > from these points 



before starting the sampling procedure. The Gillespie SSA is detailed in Table |2.1 
In the case of limit cycles, several points along the limit cycle can be used as starting 
points for the SSA. Either way, we denote the S starting positions for the trajectories 
by {yk}k=i The trajectories simulated up to some time B -\- T > using the 

SSA can help inform us about the regions of domain which will have non- negligible 
invariant density. Since the trajectories will contain many points, we cannot store 
all of the points that are simulated. Instead, we subsample from the trajectories at 
equidistant time points, at a rate Q > points per unit time. This leaves us with a 
set of sampled points 

{zMlSfJiCR^ (4.2) 
from the invariant distribution, where [QT\ denotes the integer part of the real num- 



ber QT. We hope to recover, from this set of points (4.2), information about what 
might be an optimal finite element mesh. We also denote by xf^^^ (resp. xf^^^) the 
maximal (resp. minimal) value of the i-th chemical species, i = 1,2,3, during time 
intervals [5, B -\-T] of all S simulations. These numbers will be useful in ( |4.4[ )-( [X5| ). 
In step [3], we make the approximation that the region which contains the majority 
of the invariant density, is a subset of the union of a set of ellipsoids centred at each 
point ( |4.2| ), namely 

T = [jSr{zk,i). (4.3) 

k,l 

Here, Sri^k^i) is an ellipsoid with radii r = [ri,r2,r3] centred at point z^k.i foi" ^ = 
1,2,. ..,5, / = 1,2,..., \ QT\ . These radii can be picked to be proportional to the 
range of each chemical species using the parameter /3i > and numbers x^^", x^^^, 
i = 1,2,3, computed in step [2]. Namely, we define 

grange ^ ^max _ ^min^ Z = 1, 2, 3, (4.4) 

and 

r = /3i (xf x^""^^) . (4.5) 

One should note that the parameters B and T might in principle be different for 
different initial conditions y^^ but to simplify the notation, we do not stress this fact 
by using the notation Bk and Tk and use simply B and T. 

4.2. Stage II: Automatic detection of the domain. Next, we would like 



to identify the domain on which we wish to solve (1.1). Since we already have an 
approximation of the region which contains the majority of the probability density, 
we can simply fit a cuboid around this region, and then extend it by a factor. The 
only other condition that we enforce is that the domain must be contained by the 
positive quadrant of the state space. In step [4], we pick a parameter > to define 
how much we would like to extend the cuboid: 
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(4.6) 



where 



A = (max{0, x^" - P2 ^r^'}, ^r"" + P2 X 



range \ 



1,2,3. 



We then wish to solve on the domain with boundary conditions (2.4) on dVt. 



4.3. Stage III: Construction of the adaptive mesh. In this stage, we con- 
struct the adaptively refined mesh. In the previous stage we naturahy defined the 



computational domain as a cuboid (4.6). Therefore, we base our refinement approach 



on simple splitting of a cuboid into 8 congruent sub-cuboids, as shown in Figure |4T 
However, any other standard mesh refinement technique can be utilised instead. 




Fig. 4.1. Example of a single cuboid after the proposed refinement into 8 sub-cuboids. 

An advantage of the refinement of cuboids into 8 sub-cuboids is its simplicity. A 
disadvantage is that these refinements produce so-called hanging nodes in the mesh, 
see Figure [42] (left). 

In step [5] of the mesh generation, our mesh consists of a single large cuboid. We 
then, in step [6], iteratively refine cuboids in the following fashion. In each iteration of 
the refinement procedure, we loop over all of the current cuboids that exist following 
the previous iteration of the method. For each of the cuboids Ii x I2 x /a, where 
/i, /2, /a C [0, 00) are intervals, we compute the minimum distance between all points 



in this cuboid, and the set of points (4.3) in each coordinate, namely 

dist(r,/i) = inf{|ai - bi\ : [ai,a2,as] G F, 6^ G Ii}, i = 1,2,3. 

We refine the given cuboid Ii x I2 x 1^ if 

dist(F,/^) < |/^| is satisfied for every i = 1,2,3. (4.7) 

If a cuboid is to be refined, it is split into 8 cuboids of equal volume. As detailed in 
step [7], this refinement condition can be iterated for a set number of times G N, 
where N = {1,2, 3,...}, or until the number of cuboids present in the mesh is as large 
as we would like or can deal with given our computational resources. 
Notice that the algorithm described above produces hanging nodes of order one only. 
This means that two neighbouring elements in the mesh are either of equal size or 
one of them has eight times greater volume than the other one. This is important 
for practical implementation, because the hanging nodes of order one are much easier 



to work with than the hanging nodes of higher orders. See 37 for more details and 
Figure |4.2[ b) for an illustration in 2D. In order to guarantee the continuity of the 
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Fig. 4.2. (a) 1** order hanging nodes of different types in cuboid meshes, (b) Hanging nodes 
of order 1 and order 2 on a 2D mesh. 




approximation, the value of the function at the hanging node is necessarily decided 
by the values of the function at the vertices of the less refined cube. Therefore the 
hanging nodes are not in actual fact additional degrees of freedom in the problem. 



Thus, the dimension of problem (3.4) is equal to the total number of vertices in the 
generated mesh, less the number of hanging nodes. This special treatment of hanging 
nodes actually enforces on us a mesh which is highly refined in the regions which we 
wish it to be, and then the mesh becomes gradually coarser as we move away from 
those regions. 

4.4. Tetrahedral mesh. Once the cuboid mesh has been generated by steps 



[1]— [7] of Table 4.1, we can then implement a finite element method on it. In the 
numerics shown in this paper, we further split each cuboid into 6 path tetrahedra. 
However there is nothing to say that we could not use cubic elements, but we use 
tetrahedra to simplify some implementational issues. 

We choose a refinement regime in which the elements are clustered in groups of 6 



non-obtus^ tetrahedra which together form each cuboid [t]. Figure 4.3 ^a) shows an 
example element from the refinement scheme. The idea is that if we map each cuboid 
to the unit cube, then each element is exactly of the same shape and size as all the 
others. Figure [43jb) shows how 6 of these elements can tessellate into a cube. 
A lot of refinement methods in the literature involve splitting each element into several 
smaller elements fsT. However, if this is not done in a clever way, this can lead to 
degradation of the quality of the elements. That is, some of the angles of the elements 
may become too small, leading to long thin elements, which can lead to less accuracy 
in the approximation |8j. 








Fig. 4.4. Example of a unit cube after a refinement iteration. 

In the method that we propose here, instead of splitting each element, we split each 
cuboid into 8 equally sized cuboids. Each cuboid is then split, as before, into 6 equally 
sized and shaped (after a linear transformation if not a cube) tetrahedral elements. 



Figure 4.4 shows a cube which has been refined once and been split into 48 elements. 



with 27 vertices. 

5. Parameter Selection. The algorithm has several parameters T, Q, 5, 
/32 and H whose values must be decided by the user. The parameters T, Q, 
B relate to the implementation of the Gillespie SSA, and the remaining three /32, 
H relate to the selection of the domain and the mesh, given the recorded output of 
the SSA. In this section we suggest how one might go about selecting values for these 
parameters. 

5.1. SSA parameters. The purpose of simulating the SSA trajectories, as de- 



scribed in Subsection |4.1[ is to get an indication of the regions of the domain which 
contain the majority of the invariant probability density. There are four parameters 
to consider, ^gN, T>0, (5>0 and B > 0. The parameter B represents the 
length of the simulated trajectory that we ignore at the beginning of the simulation. 
This period of simulation is simply used to ensure that the Markov chain has entered 
probabilistic equilibrium. This ensures that it is highly likely the point in state space 
where the rest of the trajectory starts has non-negligible density with respect to the 
invariant distribution. The final parameter, Q, represents the rate (per unit time) at 
which we record samples from each trajectory in the time interval [B^ B-\-T]. This pa- 
rameter is necessary since if we recorded every single state that the trajectory passed 
through, we could quickly run out of memory. 



^All six dihedral angles between its faces are less than or equal to 7r/2 
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Choosing sensible values for these parameters is not easy, since it is problem de- 
pendent. However, our aim is not to run such a long trajectory that simply using a 
histogram of the values would lead to an accurate approximation of the invariant den- 
sity, as this would completely negate the need for the rest of the algorithm. Likewise, 
if the trajectory is too short, we will have a very poor approximation of the areas of 
the domain which contain the majority of the probability density. Our aim, therefore, 
is to get as good an approximation of this region as possible with the shortest possible 
trajectories. 

As with many Monte Carlo estimates, due to the law of large numbers, the error decays 
at a rate proportional to 1/a/T. It is reasonable to expect that the approximation 
of the region which contains the majority of the invariant density should decay at a 
similar rate. Since convergence diagnostics relating to Markov chains are an open area 
with no hard and fast solution 14 , we suggest that a sample trajectory be created 
a priori to ascertain the relaxation time of the system. The parameter B can simply 
be set to be a proportion of the length T, for instance B = T/10. The parameter Q 
should be picked according to how much memory one would like to set aside to store 
the sampled points. Since we are calculating minimum distances to any point in this 
set in step [6] of the algorithm, the more points there are, the longer the algorithm 
will take to run. 

5.2. Domain and mesh generation parameters. The domain and mesh gen- 
eration parameters are given by /3i > 0, > and H gN. The parameter P2 defines 



the size of the domain as seen in (4.6). Unless the trajectories from the SSA are very 
short, the estimation of the domain should be relatively trivial, and a value of order 
1 should suffice in most instances. The value of which indicates how much error 
we estimate there to be in our estimation of the region which contains the majority 
of the invariant density, should be directly proportional to 1/Vt. Too large a value 
will lead to unnecessary refinement in some areas, while too small a value will lead 
to not enough refinement and more error. The parameter H indicates the maximum 
number of splits that the initial cubes of the mesh may undergo. In practice this can 
be found at runtime, by capping the total number of elements we allow in the mesh 
due to memory and CPU constraints. 

5.3. A numerical test of accuracy. Although the method has seven param- 
eters which have to be specified, it is relatively easy to check that the computed 
results are numerically accurate. Given the parameter set {S*, T, Q, 5, /32, i?}, we 
can compare the computed results with the results obtained with the parameter set 
{S^ 2T, 2(5, 25, 2/3i, /32, + 1}. If the result does not change significantly, then we can 
conclude that our parameters were chosen sufficiently large. Here, the parameters T, 
Q, 5, and are multiplied by 2 because increasing any of these parameters will 
increase the accuracy of the solution. In a similar way H can be increased by one, 
which means that the finest level of refinement in the mesh becomes smaller by a 
factor of two in each coordinate. However, we do not propose to modify S and f32 
for the following reasons. The number of starting points S can be determined by the 



number of stable equilibria of (4.1). In this case, the parameter S does not have to 
be changed during this a posteriori test of accuracy. Furthermore, ^2 is used only to 
choose the size of the domain, and multiplying the size of the domain by 2 would lead 



to a smaller proportion of the domain being covered by F (given by (4.3)), and as 
such a more accurate solution would not be guaranteed. If one wishes to test whether 
132 is large enough, the parameter H may have to be increased as a function of P2 in 
order to maintain the same level of refinement in F. 
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6. Implementation of the saFEM. As soon as the computational mesh is 
determined, the approximation ^ of the invariant distribution p is computed 



by the standard FEM approach 11,37. In order to compute the entries of the 
stiffness matrix ( |3.5[ ), it is necessary to use suitable quadrature routines. In numerical 
examples below, we consider chemical reactions of at most second order. Therefore, 
the diffusion D(yi) and drift v(x) coefficients are polynomials of degree at most two. 
Since the finite element space consists of piecewise linear functions, it suffices to use a 
tetrahedral quadrature rule of order three which integrates cubic polynomials exactly. 
The optimal (Gauss) quadrature rule on tetrahedra of order three has 5 points 37 . 



Since the dimension of the resulting algebraic system (3.4) can be very large, even 
for relatively coarse approximations, parallelisation of the assembly process is of 
paramount importance. Several packages exist for constructing matrices in paral- 
lel. The Portable Extensible Toolkit for Scientific Computation (PETSc) is a suite of 
data structures and routines for the scalable (parallel) solution of scientific applica- 
tions modelled by PDEs fslG] . The matrix is split into sections which are controlled by 
each processor. If the degrees of freedom are ordered in a sensible way, then Message 
Passing Interface (MPI) traffic between the processors can be kept to a minimum, 
leading to excellent scaling of processors vs. computation time |6 . 



Once the stiffness matrix (3.5) is assembled, we have to find a nontrivial solution of 
the algebraic system (3.4). Practically, we look for the eigenvector corresponding to 
the zero eigenvalue of the matrix A. To find this eigenvector in parallel, we use a 
sister package of PETSc, which is called the Scalable Library for Eigenvalue Prob- 
lem Computations (SLEPc) In particular, we used the MUltifrontal Massively 
Parallel sparse direct Solver (MUMPS) [P,^ package for the preconditioning, and a 



power method to solve the resulting eigenvalue problem 32 . 

Once we have calculated the eigenvector, it is then necessary to reconstruct the ap- 
proximated function. For this, we use the information regarding the hanging nodes 
so that we reconstruct all of the vertices on the mesh. The function in question is a 
probability density and therefore we normalise the obtained finite element solution pi^ 



such that it satisfies the condition (2.5). 



7. Numerical Results. In this section, we first use a simple example chemical 
system from 13 and implement the saFEM on a range of different mesh sizes, and 
with different algorithmic parameters. Then we present the results of saFEM for the 
Oregonator f22^. 

7.1. Convergence of the numerical method. We will study the system of 
three chemical species Xi, X2 and X3 which are subject to the following system of 
five chemical reactions [13]: 

^ Xi ^ X2 ^ X3 ^ 0. (7.1) 

Then the propensity functions are defined by 

ai(t) = kiV, a2{t) = k2Xi{t), as{t) = feXsCt), 



where V is the volume of the reactor [13 19 . We consider this system with the 



following set of non-dimensionalised parameters: 



kiV = 100, k2 = k3 = 5, k^ = k5 = 1. (7.2) 
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Fig. 7.1. (a) Convergence of the saFEM with S = I, T = 10^, Q = 0.1, B = 10^, /3i = 0.01 
and = 0.55 over a range of mesh refinements H , for the system \7.1\ . Errors are given by ( |7.3| . 
(b) Convergence of the saFEM with S = 1, Q = 10, B = 10^, /3i = 0.01, ^2 = 0.4, and H = 6 over 
a range of SSA lengths T, for the system (|7.1|> . Error is given by (|7.4| . 



As discussed in Section |5) there are seven different algorithmic parameters, whose 
values determine the accuracy and efficiency of the methods. In this section we will 
analyse the effects and convergence of the method due to altering the three most 
important of these parameters. 

7.2. Convergence of approximation. First we consider how the error in the 
approximation of the invariant distribution p decays as we refine the mesh, each time 
using the same stochastic simulation, with S = 1, T — 10^, Q = 0.1, 5 = 10^, Pi = 
0.01 and P2 = 0.55. We choose S = 1 since the corresponding ODE approximation 



given by (4.1) gives us: 

dxi 

dX2 

"dT 
dt 



100 - 5x1 + 5x2, 
5xi — 6x2 — 0.5, 

X2 - X3, 



which has a single steady state at (119.5,99.5,99.5). Therefore we pick the single 
starting point for the SSA trajectory to be at (120, 100, 100). Note that the constant 
term —0.5 in the second equation comes from the derivative of the diffusion terms as 
given by (|2.3[ ), while these derivatives sum to zero for the first and third equations. 
Figure |7.1[ a) shows the convergence of the method as the number of splits H is in- 
creased from 1 to 7, where the error is calculated by comparison with an approximation 
using a mesh which has undergone up to 8 splits: 

Ettot{H) = (^J^ |p8(x) -p,,(x)|2dx^ ^ . (7.3) 

When using a relatively coarse mesh, it is possible for the value of the approximation 
to be negative in regions. Since the function represents a probability density, which 
is strictly non- negative, this does not make sense. These negative areas can also 
complicate the normalisation of the function. We solve this problem by simply cutting 
off the negative values, i.e. we multiply the approximated function by the indicator 



function over the set where the function is non- negative, before normalising by (2.5). 
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7.3. Length of Stochastic Simulation. Next we show convergence of the 
method as T is increased, with S' = 1, Q = 10, 5 = 10^, /3i = 0.01, /^s = 0.4 
and H = S. As in the previous section, the SSA trajectory is given initial condition 
(Xi,X2,X3) = (120,100,100). 

As the system is ergodic, we would expect that as the length of stochastic simulation 
increases and the number of samples increases, that our samples becomes increasingly 
representative of the invariant density in question. This gives us more information 
about where we should be refining our mesh, which in turn should give us a more 
accurate approximation. 

Figure |7.1[ b) shows the convergence of approximation as the length of stochastic 
simulation is increased. The error as plotted in Figure 7.1 (b) is the difference 
between the approximations with varying T and the approximation calculated with 
T = 10^, namely: 



Error (T) 



/ |Pio5(x) 



■pT(x)pdx 



1/2 



(7.4) 



Note that in order to more easily compare the distributions, the domain selection for 
all of the data points was made using the longest stochastic simulation with T = 10^. 
This shows that as the simulation length is increased we see a convergence of the 
chosen mesh to one which well represents the region containing the invariant density. 

7.4. Parameter In certain situations the stochastic simulations may be ex- 
pensive and we may only be able to take a few samples from the stochastic trajectory. 
In this case we still have a large amount of uncertainty about the region which contains 
the vast majority of the invariant density. Therefore we can only hope to approximate 
this by increasing the size of the ellipsoid around each sampled point that we include 



in the region F given by (4.3). Figure 7.2 shows the convergence of the method as Pi 
is increased, with 6' = 1, T = 10^ Q = 0.1, B = 10^ (32 = 0.4 and H = 6. 




Fig. 7.2. Convergence of the saFEM with (32 = 0.4, Q = 0.1, B = T = 10^ and S 
over a range of radii of uncertainty (3i , for the system ( |7.1| > . The error is defined by ( |7.5| . 



The error as plotted in Figure 7.2 is given by: 



Error(/3i) = (^^ |pi(x) -p;3,(x)|'dx 



1/2 



(7.5) 



7.5. A numerical test of accuracy. In Section 5.3 a brief test for sufficient 
convergence was suggested, where two solutions were calculated, one with parameters 

{5,T,Q,5,/3i,/32,i/}, 
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and the second with parameters 



{S,2T,2Q,25,2/3i,/32,i/ + l}. (7.6) 

If the two solutions are close up to some tolerance, then we can be fairly sure that 
the method has converged to a reasonable degree. 



The solution for the system (7.1) calculated with the parameter set 

{S = 1, T = 10^ Q = 0.05, B = 500, /3i = 10-^ p2 = 0.45, H = 7}, 

when compared with the solution using the parameters 

{5' = 1, T = 2 X 10^ Q = 0.1, B = 10^ h = 2 x 10"^ H = 8}, 

gave error of 1.99 x 10~^, which is significantly less than the norm of the solution 
(4.51 X 10~^). This gives a relative error of 4.41 x 10~^. For ease of comparison 



the domain Q chosen using parameters (7.6) was used for both meshes. 



7.6. Oregonator. Our final example is the Oregonator |22 which is a three 
species chemical system motivated by the Belousov-Zhabotinsky reaction. It exhibits 
oscillatory behaviour and is traditionally given by the following system of ODEs: 



dxi 



kiX2 — k2XiX2 + k2,xi — 2k^x\, 



dX2 7 7 1 7 / N 

= -kiX2 - k2XiX2 + -/CcXs, (7.7) 



dx3 
dt 



2k^xi - kcXs. 



We now construct a set of reactions whose behaviour is given by (7.7) in the thermo- 
dynamic limit: 

X2 ^ Xi Xi + X2 ^ Xi ^ 2Xi + 2X3 (7.8) 



2X1 ^0 X3 > X3 > X2. 

Note that the reaction with parameter kc has been split into two reactions and 
Rq in order to get the factor of half in the equation for the rate of change of X2, 
with k^ = ke = kc/2. We consider this system with the following set of dimensionless 
parameters: 

ki = 0.3, k2 = 4000, ks = 5, k^ = 1200, kc = 0.02. (7.9) 

In this parameter regime the stochastic description of these reactions exhibits os- 
cillatory behaviour. Figure |7.3[a) shows normalised trajectories of the Oregonator 



(7.8) with parameters given by (7.9), simulated using the Gillespie SSA. This figure 
demonstrates the oscillatory behaviour of the Oregonator in this parameter regime. 
One thing of note should be mentioned at this point, regarding the ergodicity of 
this system. The zero state, at the origin, is an absorbing state for this system, 
and so trivially the invariant distribution for this system is a Dirac measure on this 
state. However, we are interested in the behaviour of this system conditioned on non- 
extinction of the species. Therefore we ensure that our domain of solution does not 
include this state, and thus the transient states involved in the oscillation behaviour 
now form the regions with non-zero invariant density. 
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Fig. 7.3. (a) Trajectories of the Oregonator ( |7.8| with parameters given by ^7.9\ , simulated 
using the Gillespie SSA. Plots show normalised trajectories, with actual molecule numbers divided 
by a time average (4:7.7, 211.5 and 2 Ax resp.). (b)-(d) Approximat ion o f the log of the marginal 
inva riant distribution (conditioned on non- extinction of species Xi) of ( |7.8| > with system parameters 
( |7.9| in the (b) xi-X2 plane, (c) xi-xs plane, (d) X2-xs plane, by the saFEM with algorithmic 
parameters given by (|7.10|. 



In the figures that fohow, these algorithmic parameters were used to approximate the 
steady state distribution (conditioned on non-extinction of species Xi): 

5 = 1, T = 10^ g = 10-^ 5 = 10^ H = 8 (T.IO) 

Since this system has a single limit cycle and is ergodic, step [1] of the algorithm (in 



Table 4.1) can be replaced by running the Gillespie SSA for a period of algorithm 
time until we are reasonably sure that the chain has entered probabilistic equilibrium, 
and using the endpoint of this simulation as the starting point for the SSA in step 
[2]. Since the system has a single stable orbit, we pick S = 1. Since a priori we do 
not know a specific point on the most likely orbit, the initial condition {xi^X2^xs) = 
(100, 100, 100) was chosen. 

Since we wish to condition on non-extinction of all of the species, we must ensure that 
Xi never drops below 1. We can do this by slightly altering the domain Q of solution 



of the FPE. We do this by replacing the first line of (4.6) by ^7 = x ^2 x ^3, where 



Al = (max{l, xf^ - p2 xT^""}. xr"" + f^2 x^^') 
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FlG. 7.4. Approxim ation of the invariant distribution (conditioned on non- extinction) of (|7.8|) 
with system parameters (|7.9|. Contours plots shown in several slices. 



The longer the simulation in the step [2] of the saFEM, the more sure we can be of 
the regions which have non-negligible probability density, and therefore the larger we 
can make ^2 and Pi. We choose Pi = 0.01 and P2 = O-l- 

Using these parameters, a mesh was created by the saFEM with 8.65 x 10^ vertices, 
out of a possible 1.70 x 10^ with up to 8 splits of each cuboid, with 1.17 x 10^ hanging 
nodes. Therefore there were a total of 7.48 x 10^ degrees of freedom, a mere 4.4% of 
the degrees of freedom that would have been required to fill the domain with cuboids 
of the finest refinement level. The mesh contained 4.79 x 10^ individual tetrahedral 
elements. 



Figure 7.4 shows several slices of the approximated invariant probability density con- 
ditioned on non-extinction of species Xi, represented by contour plots. Since visu- 
alisation of a three-dimensional function on a two-dimensional page is problematic, 
we present also in Figure 7.3 'b)-(d) the marginal densities of the approximation of 
the invariant density of the Oregonator system with parameters (7.9) in 3 different 
planes. 

To verify that the method is accurately representing the invariant distribution, we can 



test the approximation as suggested in Section |5]3) The approximation was compared 
with one with the parameters given as follows: 

5' = 1, T = 5 X 10^ Q = 5 X 10-^ B = 5x 10^ H = 7, (3i=5x 10"^ 

We omit P2 here as the same domain as was chosen for the original parameter set was 
used for ease of comparison of the two densities. The difference between the two 
solutions was 2.89 x 10~^^, with the original solution having an norm of 1.74 x 10~^, 
giving a relative error of 1.66 x 10~^, verifying that our solution is well converged. 

8. Discussion. In this paper we have presented a method for solution of the 
steady-state Fokker-Planck equation for chemical systems. The method is based on 
the use of stochastic trajectories to estimate the region in which we must refine our 
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finite element mesii tiie most. Ttiis mettiod is only valid for systems which do not 
have large regions of the domain with non-negligible invariant density that can be 
accurately approximated by a linear function. Our assumption that regions which 
have high enough probability density require high levels of refinement would fail in 
this case, and the mesh produced would be inefficient. 

The use of the Fokker-Planck equation as a means to model chemical reactions of the 
mesoscale has been present in the literature for some time 24,36 . As computational 
resources have improved, the solving of such equations in more than one dimension has 
become tractable. It has been shown that using appropriate methods (such as finite 
volume methods), the solution of the Fokker-Planck equation is far more efficient than 
using an SSA for the computation of both time-dependent and steady state solutions. 



and with a high degree of accuracy [21,36 



The Fokker-Planck equation does suffer badly from the curse of dimensionality, and 
many different methods exist which attempt to find approximations of the solution 
of these equations efficiently in more than one dimension. Some involve dimension 
reduction of the equations themselves, through simplifying assumptions, and through 



coupling with macroscopic reaction rate equations for some of the species 33 . Sparse 
grid methods have also been used in an attempt to find the invariant density of 
the chemical master equation 27 . Sparse grids attempt to overcome the curse of 
dimensionality by reducing greatly the number of degrees of freedom considered. This 
CME solution can also be coupled to Fokker-Planck equations for the more abundant 
species using a hybrid method 127 . Other methods have been used to approximate the 



solution of the chemical master equation, including adaptive wavelet compression 29 
and finite state projection fl6''34\ 

There are several other ways that stochastic trajectories could be calculated in the 
saFEM, other than the standard SSA as detailed in Table |2.1[ The r-leap method 
[26] approximates the trajectory by modelling several reactions in each iteration of 
the algorithm, up to a fixed time step. This can accelerate the simulation of the 
trajectory, but could incur some errors, including the possibility that the trajectory 
may become negative in one or more of the chemical species. Likewise, one could 
use numerical approximations of the chemical Langevin equation 25 , a stochastic 
differential equation (SDE) with corresponding Fokker-Planck equation given by (2.2), 
to create trajectories. As with the r-leap method, the trajectories from this SDE can 
become negative, and so with both of these methods proper treatment of boundary 
conditions is key to accurately approximating the region of the positive quadrant that 
contains the majority of the probability density. 

More information could also be extracted from the stochastic trajectories. A rudi- 
mentary approximation of the density could be made using the samples taken from 
the simulations, and used as a preconditioner for the eigensolver. As the eigensolver 
used is iterative, even a rough guess of the form of the solution could help to re- 
duce computation times. It has been previously shown that combining stochastic 
simulations and solutions of differential equations can be used for preconditioning of 
computations 35 . 

Acknowledgements: The research leading to these results has received funding from 
the European Research Council under the European Community's Seventh Framework 
Programme {FP7/2007-2013) / ERG grant agreement No. 239870. This publication 
was based on work supported in part by Award No KUK-C 1-0 13-04, made by King 
Abdullah University of Science and Technology (KAUST). The work of Simon Cotter 
was also partially supported by a Junior Research Fellowship of St Cross College, Uni- 



17 



versity of Oxford. Tomas Vejchodsky would like to acknowledge the support by the 
Grant Agency of the Academy of Sciences (project. IAA100760702) and the Academy 
of Sciences of the Czech Republic (institutional research plan no. AV0Z0190503). 
Radek Erban would also like to thank Somerville College, University of Oxford, for 
a Fulford Junior Research Fellowship; Brasenose College, University of Oxford, for a 
Nicholas Kurti Junior Fellowship; the Royal Society for a University Research Fellow- 
ship; and the Leverhulme Trust for a Philip Leverhulme Prize. This prize money was 
used to support a research visit of Tomas Vejchodsky in Oxford. 

REFERENCES 

[1] P. Amestoy, I. Duff, J. L'Excellent, and J. Koster. A fully asynchronous multifrontal solver 
using distributed dynamic scheduling. SIAM Journal on Matrix Analysis and Applications^ 
23(1):15"41, 2001. 

[2] P. Amestoy, A. Guermouche, J. L'Excellent, and S. Pralet. Hybrid scheduling for the parallel 

solution of linear systems. Parallel Computing^ 32(2):136-156, 2006. 
[3] A. Arkin, J. Ross, and H. McAdams. Stochastic kinetic analysis of developmental pathway 

bifurcation in phage 1-infected Escherichia coli cells. Genetics^ 149:1633-1648, 1998. 
[4] I. Babuska and W. Rheinboldt. Error estimates for adaptive finite element computations. SIAM 

Journal on Numerical Analysis^ 15(4):736-754, 1978. 
[5] S. Balay, J. Brown, K. Buschelman, W. Gropp, D. Kaushik, M. Knepley, L. Mclnnes, B. Smith, 

and H. Zhang. PETSc Web page, 2011. http://www.mcs.anl.gov/petsc. 
[6] S. Balay, W. Gropp, L. Mclnnes, and B. Smith. Efficient management of parallelism in object 

oriented numerical software libraries. In Modern Software Tools in Scientific Computing^ 

pages 163-202, 1997. 

[7] L. Beilina, S. Korotov, and M. Kfizek. Nonobtuse tetrahedral partitions that refine locally 

towards fichera-like corners. Applications of Mathematics, 50:569—581, 2005. 
[8] J. Bey. Tetrahedral grid refinement. Computing, 55:355-378, 1995. 

[9] Y. Cao, D. Gillespie, and L. Petzold. Multiscale stochastic simulation algorithm with stochastic 

partial equilibrium assumption for chemically reacting systems. Journal of Computational 

Physics, 206:395-411, 2005. 
[10] Y. Cao, D. Gillespie, and L. Petzold. The slow-scale stochastic simulation algorithm. Journal 

of Chemical Physics, 122(1):14116, 2005. 
[11] P. Ciarlet. The finite element method for elliptic problems. North-Holland Publishing Co., 

Amsterdam, 1978. Studies in Mathematics and its Applications, Vol. 4. 
[12] S. Cotter, T. Vejchodsky, and R. Erban. Stochastic bifurcations of biochemical systems, in 

preparation, 2011. 

[13] S. Cotter, K. Zygalakis, I. Kevrekidis, and R. Erban. A constrained approach to multi- 
scale stochastic simulation of chemically reacting systems. Journal of Chemical Physics, 
135:094102, 2011. 

[14] M. Cowles and B. Carlin. Markov chain monte carlo convergence diagnostics: A comparative 
review. Journal of the American Statistical Association, 91(434):883-904, 1996. 

[15] E. Doedel, A. Champneys, T. Fairgrieve, Y. Kuznetsov, B. Sandstede, and X. Wang. AUTO 97: 
Continuation and bifurcation software for ordinary differential equations (with homcont), 
1997. 

[16] B. Drawert, M. Lawson, L. Petzold, and M. Khammash. The diffusive finite state projec- 
tion algorithm for efficient simulation of the stochastic reaction-diffusion master equation. 
Journal of Chemical Physics, 132:074101, 2010. 

[17] W. E, D. Liu, and E. Vanden-Eijnden. Nested stochastic simulation algorithm for chemical 
kinetic systems with disparate rates. Journal of Chemical Physics, 123:194107, 2005. 

[18] R. Erban, S. J. Chapman, I. Kevrekidis, and T. Vejchodsky. Analysis of a stochastic chemical 
system close to a sniper bifurcation of its mean-field model. SIAM Journal on Applied 
Mathematics, 70(3):984-1016, 2009. 

[19] R. Erban, S. J. Chapman, and P. Maini. A practical guide to stochastic simulations of reaction- 
diffusion processes. 35 pages, available as http://arxiv.org/abs/0704.1908, 2007. 

[20] R. Erban, I. Kevrekidis, D. Adalsteinsson, and T. Elston. Gene regulatory networks: A coarse- 
grained, equation-free approach to multiscale computation. Journal of Chemical Physics, 
124:084106, 2006. 

[21] L. Perm, P. Lotstedt, and P. Sjoberg. Conservative solution of the fokker-planck equation for 



18 



stochastic chemical reactions. BIT Numerical Analysis, 46:61-83, 2006. 
[22] R. Field and R. Noyes. Oscillations in chemical systems. IV. limit cycle behavior in a model of 

a real chemical reaction. Journal of Chemical Physics, 60(5): 1877-1884, 1974. 
[23] D. Gillespie. Exact stochastic simulation of coupled chemical reactions. Journal of Physical 

Chemistry, 81(25):2340-2361, 1977. 
[24] D. Gillespie. Approximating the master equation by Fokker-Planck equations for single-variable 

chemical systems. Journal of Chemical Physics, 72(10) :5363-5370, 1980. 
[25] D. Gillespie. The chemical Langevin equation. Journal of Chemical Physics, 113(l):297-306, 

2000. 

[26] D. Gillespie. Approximate accelerated stochastic simulation of chemically reacting systems. 

Journal of Chemical Physics, 115(4):1716-1733, 2001. 
[27] M. Hegland, A. Hellander, and P. Lotstedt. Sparse grids and hybrid methods for the chemical 

master equation. BIT Numerical Mathematics, 48(2):265-283, 2008. 
[28] V. Hernandez, J. Roman, and V. Vidal. SLEPc: A scalable and flexible toolkit for the solution 

of eigenvalue problems. ACM Transactions on Mathematical Software, 31(3):351-362, 

2005. 

[29] T. Jahnke and T. Udrescu. Solving chemical master equations by adaptive wavelet compression. 

Journal of Computational Physics, 229(16):5724 - 5741, 2010. 
[30] S. Kar, W. Baumann, M. Paul, and J. Tyson. Exploring the roles of noise in the eukaryotic 

cell cycle. Proceedings of the National Academy of Sciences USA, 106:6471-6476, 2009. 
[31] S. Korotov and M. Kek. Acute type reflnements of tetrahedral partitions of polyhedral domains. 

SIAM Journal on Numerical Analysis, 39(2):724-733, 2002. 
[32] V. Kublanovskaya. On some algorithms for the solution of the complete eigenvalue problem. 

U.S.S.R. Computational Mathematics and Mathematical Physics, l(3):555-570, 1961. 
[33] P. Lotstedt and L. Perm. Dimensional reduction of the fokker-planck equation for stochastic 

chemical reactions. Multiscale Modelling and Simulation, 5:593-614, 2006. 
[34] B. Munsky and M. Khammash. The flnite state projection algorithm for the solution of the 

chemical master equation. The Journal of Chemical Physics, 124(4) :044104, 2006. 
[35] L. Qiao, R. Erban, C. Kelley, and I. Kevrekidis. Spatially distributed stochastic systems: 

Equation-free and equation-assisted preconditioned computation. Journal of Chemical 

Physics, 125:204108, 2006. 
[36] P. Sjoberg, P. Lotstedt, and J. Elf. Fokker-planck approximation of the master equation in 

molecular biology. Computing and Visualization in Science, 12:37-50, 2009. 
[37] P. Solm, J. Cerveny, and I. Dolezel. Arbitrary-level hanging nodes and automatic adaptivity 

in the ^p-FEM. Mathematics and Computers in Simulation, 77(1):117-132, 2008. 
[38] R. Verfiirth. A review of a posteriori error estimation and adaptive mesh-refinement techniques. 

Ghichester: John Wiley & Sons. Stuttgart: B. G. Teubner., 1996. 
[39] J. Villar, H. Kueh, N. Barkai, and S. Leibler. Mechanisms of noise-resistance in genetic oscil- 
lators. Proceedings of the National Academy of Sciences USA, 99(9):5988-5992, 2002. 
[40] C. Zenger. Sparse grids. In W. Hackbusch, editor. Parallel algorithms for partial differential 

equations, volume 31 of Notes on Numerical Fluid Mechanics, pages 241-251. Vieweg- 

Verlag, 1991. 



19 



